#!/bin/bash -e

function info() {
echo Usage: `basename $0` '[-s name] [-t threads] [-2] [-b build] [-g ref_genome] [-r [repos]] reads1.fq reads2.fq'
exit 65
}

while getopts  ":b:s:t:p:2g:r" opts
do
        case  $opts  in
		b) asemmbly=$OPTARG;;
        s) sample_name=$OPTARG;;
		p) out_prefix=$OPTARG;;
		t) threads=$OPTARG;;
		g) ref_genome=$OPTARG;;
		r) repos=T;;
		2) tool=bt;;
		\?) echo info;;
        esac
done
shift $(($OPTIND - 1))


if [ $# -lt 2 ]; then info; fi

if test -z $threads; then threads=2; fi

. $var

# clipping_penalty=1000
clipping_penalty=5


echo pair-end reads 1: $1
echo pair-end reads 2: $2
echo sample: $sample_name. threads: $threads.

if test "$tool" = "bt"
then

echo;echo;echo bowtie2 alignment
bowtie2 \
-x $ref_genome \
-1 $1 \
-2 $2 \
-S $out_prefix.sam \
--local \
--rg-id $sample_name \
--rg SM:$sample_name \
--rg PU:$sample_name \
--rg LB:$sample_name \
--rg PL:ILLUMINA \
-p $threads

else

echo;echo;echo start bwa mem ...
$bwa mem -M -L $clipping_penalty -t$threads -R $read_group $ref_genome \
$1 \
$2 \
> $out_prefix.sam

fi

if test -n "$repos" -o "$ref_genome" != "$data_path/ref/$genome_assembly/$genome_assembly.fa"; then
samtools view $out_prefix.sam|perl -ne 'my @fs=split "\t"; if ($fs[2] != "*") {my @pos= split ":", $fs[2]; $fs[2]=$pos[0]; $fs[3]+=$pos[1]; $fs[7]+=$pos[1]; $"="\t"; print "@fs"}' >$out_prefix.repos.sam

samtools view -H $out_prefix.sam|grep '@RG'|cat $data_path/ref/$genome_assembly/sam.header.txt -  $out_prefix.repos.sam > $out_prefix.repos.header.sam && rm $out_prefix.repos.sam
fi

. $cmd_done

